home *** CD-ROM | disk | FTP | other *** search
/ MacWorld 1997 September / Macworld (1997-09).dmg / Serious Software / Cherwell Scientific Demos / pro Fit / pro Fit 5.0 demo (fpu).sea / pro Fit 5.0 demo (fpu) / Examples / Plotting - vector fields / DipoleField program < prev    next >
Text File  |  1996-04-15  |  3KB  |  120 lines
  1. {
  2. This program calculates the electric field of two opposite
  3. charges arranged at points -1/0 and 1/0.
  4. To run this program, first add it to pro Fit by choosing
  5. "Add to Menu" from the "Misc" menu. Then open a new
  6. drawing window. Then select the program ("DipoleField")
  7. from the "Run Program" submenu in the "Calc" menu.
  8.  
  9. This program can be adapted to plotting other vector fields by
  10. modifying the procedure "GetField" only.
  11. }
  12.  
  13. program DipoleField;
  14.  
  15. const 
  16.       arrowAngle = 30; { arrow opening angle }
  17.       arrowLen = 0.02; { arrow length }
  18.       rangeX = 2;
  19.       rangeY = 1;
  20.       zoom = 0.02;     { zoom factor for drawing vectors }
  21.       maxArrow = 0.10; { the max length of an arrow }
  22.  
  23. var  Ex, Ey: real;     { the results of ElectricField }
  24.      E1x, E1y: real;
  25.      
  26.      x,y: real;
  27.      VectorX, VectorY:real; {set by the function GetField }
  28.      sinA, cosA: real; { sine and cos of arrow angle times arrow length }
  29.  
  30.  
  31.  
  32. procedure DrawArrow(x,y, vx,vy);
  33.   { draws an arrow at position x,y with vector vx, vy }
  34.   var tipX, tipY: real;
  35.       ax, ay, r: real;
  36. begin
  37.   vx := vx*zoom; vy := vy*zoom;
  38.   r := sqrt(sqr(vx) + sqr(vy));
  39.   if (r < maxArrow) and (r > 1e-30) then
  40.   begin
  41.    tipX := x + vx/2;
  42.    tipY := y + vy/2;
  43.    MoveTo(x-vx/2, y-vy/2);
  44.    LineTo(tipX, tipY);
  45.    if r >= arrowLen then   { show arrow only if shaft long enough }
  46.    begin
  47.     vx := vx/r; vy := vy/r;
  48.     ax := cosA*vx - sinA*vy;
  49.     ay := sinA*vx + cosA*vy;
  50.     MoveTo(tipX-ax, tipY-ay);
  51.     LineTo(tipX, tipY);
  52.     ax := cosA*vx + sinA*vy;
  53.     ay := -sinA*vx + cosA*vy;
  54.     LineTo(tipX-ax, tipY-ay);
  55.    end;
  56.   end;
  57. end; { DrawArrow }
  58.  
  59.  
  60. procedure ElectricField(x,y, xx0, yy0, c);
  61.   { returns the field of a single charge }
  62.  var r3;
  63. begin
  64.  x := x-xx0;
  65.  y := y-yy0;
  66.  r3 := (sqr(x) + sqr(y)) ^1.5;
  67.  if (r3 < 1e-30) then
  68.  begin
  69.    Ex := 0; Ey := 0;
  70.  end
  71.  else
  72.  begin
  73.    Ex := c * x/r3; Ey := c * y/r3;
  74.  end;
  75. end;
  76.  
  77. {*************************************************************************}
  78. { this function calculates the vector field to be plotted. }
  79. { substitute it with a function of yours if you want to plot }
  80. { other vector rields. Presently it calls the procedure     }
  81. { ElectricField to calculate the field of a dipole.        }
  82.  
  83. procedure GetField(x,y:real);
  84. begin
  85.       ElectricField(x,y, 1, 0, 1);
  86.       E1x := Ex; E1y := Ey;
  87.       ElectricField(x,y, -1, 0, -1);
  88.       
  89.                         VectorX:= Ex+E1x;    
  90.                         VectorY:= Ey+E1y
  91. end;
  92. {*************************************************************************}
  93.  
  94.  
  95. begin
  96.   sinA := sin(arrowAngle*π/180) * arrowLen;
  97.   cosA := cos(arrowAngle*π/180) * arrowLen;
  98.   SetNewGraphRect(30, 30, 430, 230);
  99.   CreateNewGraf(-rangeX,rangeX,-rangeY,rangeY,0,0);
  100.   SetLineStyle(0.5, 1);
  101.   OpenCurve('electric field');
  102.   x := -rangeX;
  103.   while x <= rangeX do
  104.   begin
  105.     y := -rangeY;
  106.     while y <= rangeY do
  107.     begin
  108.       GetField(x,y); {This function calculates the vector from the}
  109.                      {coordinate x,y and sets VectorX,VectorY}
  110.       DrawArrow(x,y, VectorX,VectorY);
  111.       y := y+rangeY/10;
  112.     end; {while}
  113.     x := x+rangeY/10;
  114.   end; {while}
  115.   CloseCurve;
  116.   SetLineStyle(1, 1); {reset line style}
  117. end;
  118.  
  119.   
  120.